High dimensional single-cell analysis reveals iNKT cell developmental trajectories and effector fate decision¶
In [1]:
import numpy as np
import pandas as pd
import scanpy as sc
import mnnpy as mnnpy
import numba as numba
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
import logging
mpl_logger = logging.getLogger('matplotlib')
mpl_logger.setLevel(logging.WARNING)
sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_versions()
In [2]:
sc.settings.set_figure_params(dpi=130)
Loading 10x data¶
In [3]:
adata = sc.read_h5ad('./output/adata.preprocessing.h5ad')
adata.raw = adata
print("adata shape = ", adata.shape)
adata.obs['sample']
adata.obs['batch'] = adata.obs['sample']
adata.obs
Out[3]:
In [4]:
sc.pl.scatter(adata, x='n_counts', y='percent_mito', color="sample")
sc.pl.scatter(adata, x='n_counts', y='percent_ribo', color="sample")
sc.pl.scatter(adata, x='n_counts', y='n_genes', color="sample")
In [5]:
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.25, batch_key='sample')
sc.pl.highly_variable_genes(adata)
In [6]:
adata_wt = adata[adata.obs['sample'].isin(['wt'])]
adata_ko = adata[adata.obs['sample'].isin(['ko'])]
print("wt shape = ", adata_wt.shape)
print("ko shape = ", adata_ko.shape)
In [7]:
hvgs = adata.var.loc[adata.var['highly_variable_intersection'] == True, :].index
hvgs
Out[7]:
In [8]:
corrected = mnnpy.mnn_correct(adata_wt, adata_ko, var_subset=hvgs, save_raw= True, batch_categories = ['wt','ko'])
In [9]:
adata = corrected[0]
adata.write('./output/adata.mnncorrect.new.h5ad')
In [10]:
def get_cluster_proportions(adata,
cluster_key="cluster_final",
sample_key="replicate",
drop_values=None):
"""
Input
=====
adata : AnnData object
cluster_key : key of `adata.obs` storing cluster info
sample_key : key of `adata.obs` storing sample/replicate info
drop_values : list/iterable of possible values of `sample_key` that you don't want
Returns
=======
pd.DataFrame with samples as the index and clusters as the columns and 0-100 floats
as values
"""
adata_tmp = adata.copy()
sizes = adata_tmp.obs.groupby([cluster_key, sample_key]).size()
props = sizes.groupby(level=1).apply(lambda x: 100 * x / x.sum()).reset_index()
props = props.pivot(columns=sample_key, index=cluster_key).T
props.index = props.index.droplevel(0)
props.fillna(0, inplace=True)
if drop_values is not None:
for drop_value in drop_values:
props.drop(drop_value, axis=0, inplace=True)
return props
def plot_cluster_proportions(cluster_props,
cluster_palette=None,
xlabel_rotation=0):
fig, ax = plt.subplots(dpi=130)
fig.patch.set_facecolor("white")
cmap = None
if cluster_palette is not None:
cmap = sns.palettes.blend_palette(
cluster_palette,
n_colors=len(cluster_palette),
as_cmap=True)
cluster_props.plot(
kind="bar",
stacked=True,
ax=ax,
legend=None,
colormap=cmap
)
ax.legend(bbox_to_anchor=(1.01, 1), frameon=False, title="Cluster")
sns.despine(fig, ax)
ax.tick_params(axis="x", rotation=xlabel_rotation)
ax.set_xlabel(cluster_props.index.name.capitalize())
ax.set_ylabel("Proportion")
fig.tight_layout()
return fig
In [11]:
#adata = sc.read_h5ad('./output/adata.mnncorrect.h5ad')
print("adata shape = ", adata.shape)
adata.obs
Out[11]:
In [12]:
adata.var['highly_variable'] = adata.var['highly_variable_intersection']
In [13]:
sc.tl.pca(adata, svd_solver='arpack')
In [14]:
sc.pl.pca_variance_ratio(adata, log = False)
In [15]:
sc.pp.neighbors(adata)
In [16]:
sc.tl.umap(adata)
In [17]:
sc.pl.umap(adata, color=['sample','batch'], edges = False)
In [18]:
#metadata = pd.read_csv('preprocessing_metadata.tsv', sep = '\t')
#adata.obs['cell_type'] = metadata['cell_type']
sc.tl.louvain(adata, resolution=1.5)
sc.pl.umap(adata, color=['louvain'])
In [19]:
prop = get_cluster_proportions(adata,"sample","louvain")
prop
Out[19]:
In [20]:
plot = plot_cluster_proportions(prop)
In [21]:
sc.tl.rank_genes_groups(adata, 'louvain')
In [22]:
pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(10)
Out[22]:
In [23]:
sc.pl.umap(adata, color=['Mki67', 'Top2a', 'Ube2c'], color_map="jet") # Cycling
In [24]:
sc.pl.umap(adata, color=['Lef1','Itm2a','Ccr9','Id3','Ldhb'], color_map="jet") # NKT0
In [25]:
sc.pl.umap(adata, color=['Zbtb16','Plac8','Il4','Icos','Izumo1r'], color_map="jet") # NKT2
In [26]:
sc.pl.umap(adata, color=['Fhl2','Ly6a'], color_map="jet") # NKT1a
In [27]:
sc.pl.umap(adata, color=['Il1r1', 'Rorc','Ccr6','Emb'], color_map="jet") # NKT17
In [28]:
sc.pl.umap(adata, color=['Id2','Nkg7','Klrb1c','Klrd1','Il2rb','Fcer1g'], color_map="jet") # NKT1
In [29]:
sc.pl.umap(adata, color=['Ifit1','Ifit3','Isg15'], color_map="jet") # NKT1d
In [30]:
sc.settings.set_figure_params(dpi=130)
sc.pl.umap(adata, color=['louvain'], legend_loc = 'on data', legend_fontsize = 6)
In [31]:
marker_genes = ['Lef1','Itm2a','Mki67','Top2a','Zbtb16','Plac8','Izumo1r','Ifit1','Ifit3','Isg15','Ly6a','Fhl2','Nkg7','Fcer1g','Gzma','Gzmb','Rorc','Serpinb1a','Tmem176a','Tmem176b']
ax = sc.pl.dotplot(adata, marker_genes, groupby='louvain')
In [32]:
mean_cellType = np.empty((len(adata.obs['louvain'].unique()), adata.raw.shape[1]),
dtype=float, order='C')
In [33]:
raw_adata = adata.raw.X.toarray()
In [34]:
for i in range(0, len(adata.obs['louvain'].unique())):
#print(adata.obs['phenograph'].unique()[i])
mean_cellType[i,:] = np.mean(raw_adata[adata.obs['louvain'] == adata.obs['louvain'].unique()[i], :], axis = 0)
In [35]:
mean_df = pd.DataFrame(np.corrcoef(mean_cellType), index = adata.obs['louvain'].unique(), columns = adata.obs['louvain'].unique())
In [36]:
import seaborn
ax = seaborn.clustermap(mean_df, cmap="jet").fig.suptitle('louvain')
In [41]:
new_cluster_names = {
'0':'NKT1b', '1':'NKT1b', '2':'NKT1b', '3':'NKT2b', '4':'NKT1c',
'5':'NKT1c','6':'NKT2b', '7':'NKT17', '8':'NKT2a', '9':'NKT1b', '10':'NKT1d',
'11':'Cycling NKT','12':'NKT1a', '13':'NKT0'}
vect = []
for i in range(0, len(adata.obs['louvain'])):
vect = vect + [new_cluster_names[str(adata.obs['louvain'][i])]]
adata.obs['cell_type'] = vect
adata.uns['cell_type_colors'] = [
'#F3766E', # Cycling NKT
'#A31E22', # NKT0
'#50C878', # NKT1a
'#2E8B57', # NKT1b
'#0B6623', # NKT1c
'#4CBB17', # NKT1d
'#FEC85A', # NKT2a
'#FAE600', # NKT2b
'#2da9d2'] # NKT17
In [42]:
sc.pl.umap(adata, color='cell_type', title='', frameon=False)
In [43]:
adata.obs['cell_type'].cat.categories
adata.obs['cell_type'].cat.reorder_categories(['NKT0','Cycling NKT','NKT17','NKT2a','NKT2b','NKT1a','NKT1b','NKT1c','NKT1d'], inplace = True)
adata.uns['cell_type_colors'] = [ '#A31E22', # NKT0
'#F3766E', # Cycling NKT
'#2da9d2', # NKT17
'#FEC85A', # NKT2a
'#FAE600', # NKT2b
'#50C878', # NKT1a
'#2E8B57', # NKT1b
'#0B6623', # NKT1c
'#4CBB17'] # NKT1d
sc.pl.umap(adata, color='cell_type', title='', frameon=False, save='figure.8e.pdf')
In [44]:
prop = get_cluster_proportions(adata,"sample","cell_type")
prop
Out[44]:
In [45]:
plot = plot_cluster_proportions(prop, cluster_palette=None, xlabel_rotation=90)
In [51]:
prop = get_cluster_proportions(adata,"cell_type","sample")
prop
Out[51]:
In [48]:
colors = ['#A31E22','#F3766E','#2da9d2','#FEC85A','#FAE600','#50C878','#2E8B57','#0B6623','#4CBB17']
plot = plot_cluster_proportions(prop, cluster_palette=colors, xlabel_rotation=90)
plot.savefig("figure.7e.pdf", dpi=300)
In [52]:
marker_genes = ['Lef1','Itm2a','Mki67','Top2a','Zbtb16','Plac8','Izumo1r','Ifit1','Ifit3','Isg15','Ly6a','Fhl2','Nkg7','Fcer1g','Gzma','Gzmb','Rorc','Serpinb1a','Tmem176a','Tmem176b']
ax = sc.pl.dotplot(adata, marker_genes, groupby='cell_type')
In [53]:
adata.obs['umap_1'] = pd.DataFrame(adata.obsm['X_umap']).iloc[:,0].tolist()
adata.obs['umap_2'] = pd.DataFrame(adata.obsm['X_umap']).iloc[:,1].tolist()
adata.obs.to_csv(path_or_buf = 'output/metadata.aggr.tsv', sep = '\t', index = True)
In [54]:
adata.obsm
Out[54]:
In [55]:
adata.write('./output/aggr.ann.h5ad')
In [56]:
sc.external.exporting.cellbrowser(adata,
'cellbrowser/aggr',
'aggr',
annot_keys=['cell_type', 'louvain', 'sample',
'n_counts', 'n_genes', 'percent_mito', 'percent_ribo'],
cluster_field='cell_type', nb_marker=30,
skip_matrix=False, do_debug=True)
In [57]:
color_dictionnary = {
'Cycling NKT':'#F3766E',
'NKT0':'#A31E22',
'NKT1a':'#50C878',
'NKT1b':'#2E8B57',
'NKT1c':'#0B6623',
'NKT1d':'#4CBB17',
'NKT2a':'#FEC85A',
'NKT2b':'#FAE600',
'NKT17': '#2da9d2',
}
df = pd.DataFrame.from_dict(color_dictionnary, orient = 'index', columns = ['color'])
df.to_csv(path_or_buf = 'cellbrowser/aggr/clusterColor.tsv', sep = '\t', index_label = 'Cluster')
density plot¶
In [58]:
sc.tl.embedding_density(adata, basis='umap', groupby='sample')
In [59]:
adata.obs['sample'].value_counts()
sc.pl.embedding_density(adata, basis='umap', key='umap_density_sample', group='wt', save="wt.pdf")
sc.pl.embedding_density(adata, basis='umap', key='umap_density_sample', group='ko', save="ko.pdf")
